
#ifndef UTILSCUH
#define UTILSCUH
// ============= Global configurations ==============
#include <cuComplex.h>
#include <complex>
#include <iomanip>
#include <algorithm>

#define BLOCKSIZE 1024
using Real = float;
using Comp = ::std::complex<Real>;
using namespace std;


// constexpr Real Pi = 3.1415926535897932384626433832795;
// constexpr Real _2Pi = 6.283185307179586476925286766559;

// constexpr Real Pi = 3.14f;
// constexpr Real _2Pi = 6.28f;

constexpr Real Pi = 3.1415927f;
constexpr Real _2Pi = 6.2831854f;
// =============== Windows version ctz =================
#ifdef _MSC_VER
#include <intrin.h>
typedef unsigned long DWORD;
uint32_t __inline __builtin_ctzll(uint32_t value ){
    DWORD trailing_zero = 0;

    if ( _BitScanForward( &trailing_zero, value ) ) // 从value的最低有效位向最高有效位搜索，给出1所在的index
    {
        return trailing_zero;
    }
    else
    {
        // This is undefined, I better choose 32 than 0
        return 32;
    }
}
#endif
// =====================================================

/*
    handle the string and char
*/
vector<string> split(const string& str, const string& delim) {  
	vector<string> res;  
	if("" == str) return res;  
	//先将要切割的字符串从string类型转换为char*类型  
	char * strs = new char[str.length() + 1] ; //不要忘了  
	strcpy(strs, str.c_str());   
 
	char * d = new char[delim.length() + 1];  
	strcpy(d, delim.c_str());  
 
	char *p = strtok(strs, d);  
	while(p) {  
		string s = p; //分割得到的字符串转换为string类型  
		res.push_back(s); //存入结果数组  
		p = strtok(NULL, d);  
	}  
 
	return res;  
}


template<class I>
 I ceil2(I x) {
    --x;
    x |= x >> 1;
    x |= x >> 2;
    x |= x >> 4;
    if (sizeof(I) > 1)
        x |= x >> 8;
    if (sizeof(I) > 2)
        x |= x >> 16;
    if (sizeof(I) > 4)
        x |= x >> 32;
    return x + 1;
}

template<class I>
constexpr bool isPow2(I x) { return ((x ^ (x - 1)) >> 1) == x - 1; }

Real error( Comp* A,  Comp* B, std::size_t N) {
    using namespace std;
    Real res = 0;
    for (size_t i = 0; i < N; ++i){

        res = max(res, abs(A[i] - B[i]));
        // cout << res << endl;
    }
    return res;
}

class MyException : public exception
{
public:
    MyException() : exception() {}       
};

Real error( Real* A,  Real* B, std::size_t N) {
    using namespace std;
    Real res = 0;
    Real tmp = 0;
    for (size_t i = 0; i < N; ++i){
        tmp = res;
        res = max(res, abs(A[i] - B[i]));
        
         if (tmp!=res){
         cout << fixed << setprecision(10) << A[i] << endl;
         cout << fixed << setprecision(10) << B[i] << endl;
         }
        
    }
    return res;
}

Real modulus2(const Comp A, const Comp B) {
    using namespace std;
    return abs(A - B);
}


std::mt19937_64 R(std::random_device{}());

Real randInt(int lo, int up) {
    return std::uniform_int_distribution<int>(lo, up)(R);
}

Real randReal(Real lo, Real up) {
    return std::uniform_real_distribution<Real>(lo, up)(R);
}

Comp randComp(Real lo, Real up) {
    return {randReal(lo, up), randReal(lo, up)};
}

template <typename T> 
std::complex<T> const multiply(Comp & a, Comp & b){
    return {a.real() * b.real() - a.imag() * b.imag(),
        a.imag() * b.real() + a.real() * b.imag()};
}


/* 后继 */
int getNext(int i, int m, int n)
{
    return (i%n)*m + i/n;
}

/* 前驱 */
int getPre(int i, int m, int n)
{
    return (i%m)*n + i/m;
}

/* 处理以下标i为起点的环 */
void movedata(Comp* mtx, int i, int m, int n)
{
    Comp temp = mtx[i];  // 暂存
    int cur = i;       // 当前下标
    int pre = getPre(cur, m, n);
    while(pre != i)
    {
        mtx[cur] = mtx[pre];
        cur = pre;
        pre = getPre(cur, m, n);
    }
    mtx[cur] = temp;
}

/* 转置，即循环处理所有环 */
void transpose(Comp* mtx, int m, int n)
{
    for(int i=0; i<m*n; ++i)
    {
        int next = getNext(i, m, n);
        while(next > i) // 若存在后继小于i说明重复
            next = getNext(next, m, n);
        if(next == i)   // 处理当前环 
            movedata(mtx, i, m, n);
    }
}


void transReal2Comp(Real* in,Comp* out, size_t L){
    for (int i = 0 ; i < L ; i++){

        out[i].real(in[i]);
        out[i].imag(0);
    }
}
void transComp2Real(Comp* in , Real* out,size_t L){
    for (int i = 0 ; i < L ; i++){
        out[i] = in[i].real();
    }
}

void transReal2Comp(Real* in,Comp* out, size_t h, size_t w, size_t ori_h , size_t ori_w){
       for (int i = 0 ; i < ori_h ; i++){
            for (int j = 0 ; j < ori_w ; j++){
                    out[i * w + j].real(in[i * ori_w + j]);
                    out[i * w + j].imag(0.0);
            }
        }
}
void transReal2CompWithMemset(Real* in,Comp* out, size_t n , size_t c, size_t h, size_t w, size_t ori_h , size_t ori_w){
    for (int k = 0 ; k < n ; k++ ){  
        for ( int l = 0 ; l < c ; l++ ){  
            for (int i = 0 ; i < h ; i++){
	        for (int j = 0 ; j < w ; j++){
	            if (i < ori_h && j < ori_w){
		          out[k * c * h * w + l * h * w + i * w + j].real(in[k * c * ori_h * ori_w + l * ori_h * ori_w + i * ori_w + j]);
		          out[k * c * h * w + l * h * w + i * w + j].imag(0.0);
	            } else {
	                  out[k * c * h * w + l * h * w + i * w + j].real(0.0);
		          out[k * c * h * w + l * h * w + i * w + j].imag(0.0);
						        }
		    }
		 }
	    }
        }
     }
void transReal2CompWithMemset(Real* in,Comp* out, size_t n , size_t c, size_t d, size_t h, size_t w, size_t ori_d ,size_t ori_h , size_t ori_w){
    for (int k = 0 ; k < n ; k++ ){  
        for ( int l = 0 ; l < c ; l++ ){  
            for ( int m = 0 ; m < d ; m++ ){  
                for (int i = 0 ; i < h ; i++){
	            for (int j = 0 ; j < w ; j++){
	                if (i < ori_h && j < ori_w && m < ori_d){
		            out[k * c * d * h * w + l * d * h * w + m * h * w + i * w + j].real(in[k * c * ori_d * ori_h * ori_w + l * ori_d * ori_h * ori_w + m * ori_h * ori_w + i * ori_w + j]);
		            out[k * c * d * h * w + l * d * h * w + m * h * w + i * w + j].imag(0.0);
	                } else {
	                    out[k * c * d * h * w + l * d * h * w + m * h * w + i * w + j].real(0.0);
		            out[k * c * d * h * w + l * d * h * w + m * h * w + i * w + j].imag(0.0);
						        }
		    }
		}
	    }
        }
     }
}
void transReal2CompWithMemset(Comp* in,Comp* out, size_t n , size_t c, size_t d, size_t h, size_t w, size_t ori_d ,size_t ori_h , size_t ori_w){
    for (int k = 0 ; k < n ; k++ ){  
        for ( int l = 0 ; l < c ; l++ ){  
            for ( int m = 0 ; m < d ; m++ ){  
                for (int i = 0 ; i < h ; i++){
	            for (int j = 0 ; j < w ; j++){
	                if (i < ori_h && j < ori_w && m < ori_d){
		            out[k * c * d * h * w + l * d * h * w + m * h * w + i * w + j].real(in[k * c * ori_d * ori_h * ori_w + l * ori_d * ori_h * ori_w + m * ori_h * ori_w + i * ori_w + j].real());
		            out[k * c * d * h * w + l * d * h * w + m * h * w + i * w + j].imag(in[k * c * ori_d * ori_h * ori_w + l * ori_d * ori_h * ori_w + m * ori_h * ori_w + i * ori_w + j].imag());
	                } else {
	                    out[k * c * d * h * w + l * d * h * w + m * h * w + i * w + j].real(0.0);
		            out[k * c * d * h * w + l * d * h * w + m * h * w + i * w + j].imag(0.0);
						        }
		    }
		}
	    }
        }
     }
}
void transReal2CompWithMemset(Real* in,Comp* out, size_t h, size_t w, size_t ori_h , size_t ori_w){
    for (int i = 0 ; i < ori_h ; i++){
            for (int j = 0 ; j < ori_w ; j++){
                    if (i < ori_h && j < ori_w){
		      out[i * w + j].real(in[i * ori_w + j]);
                      out[i * w + j].imag(0.0);
		    } else {
		      out[i * w + j].real(0.0);
                      out[i * w + j].imag(0.0);
		    }
            }
        }

}

void transReal2Comp(Real* in,Comp* out, size_t d, size_t h, size_t w, size_t ori_d , size_t ori_h , size_t ori_w){
    for (int k = 0 ; k < d ; k++){
        if(k < ori_d){
            for (int i = 0 ; i < h ; i++){
                if (i < ori_h){
                    for (int j = 0 ; j < w ; j++){
                        if (j < ori_w) {
                            out[k * h * w + i * w + j].real(in[k * ori_h * ori_w + i * ori_w + j]);
                            out[k * h * w + i * w + j].imag(0.0);
                        }
                    }
                }
            }
        }
    }
}

void matrixTrimPadding(Comp* in, Comp* out, size_t batch, size_t output_channel, size_t in_h, size_t in_w, size_t out_h, size_t out_w){
    if (in_h < out_h || in_w < out_w 
     //   || (in_h - out_h) % 2 != 0 
    // || (in_w - out_w) % 2 != 0
){
        std::cout << "param error" << std::endl;
        throw "Error params in matrixTrimPadding! Please check!";
        
    }
    int start_row_index = (in_h - out_h) / 2;
    int start_col_index = (in_w - out_w) / 2;
    
    for (int n = 0 ; n < batch ; n++){
        for (int k = 0 ; k < output_channel ; k++){
            for (int i = 0 ; i < out_h ; i ++){
                for (int j = 0 ; j < out_w ; j ++){
                    out[n * output_channel * out_h * out_w + k * out_h * out_w + i * out_w + j] 
                    = in[n * output_channel * in_h * in_w  + k * in_h * in_w + + (i + start_row_index) * in_w + j + start_col_index];
                }
            }
        }
    }
    
}
void matrixTrimPadding(Comp* in, Comp* out, size_t batch, size_t output_channel, size_t in_h, size_t in_w, size_t out_h, size_t out_w, int tmp_out_h_ori, int tmp_out_w_ori, int nouse){
    if (in_h < out_h || in_w < out_w 
     //   || (in_h - out_h) % 2 != 0 
    // || (in_w - out_w) % 2 != 0
){
        std::cout << "param error" << std::endl;
          // throw "Error params in matrixTrimPadding! Please check!";
        
    }
    int start_row_index = (tmp_out_h_ori - out_h) / 2;
    int start_col_index = (tmp_out_w_ori - out_w) / 2;
    
    for (int n = 0 ; n < batch ; n++){
        for (int k = 0 ; k < output_channel ; k++){
            for (int i = 0 ; i < out_h ; i ++){
                for (int j = 0 ; j < out_w ; j ++){
                    out[n * output_channel * out_h * out_w + k * out_h * out_w + i * out_w + j] 
                    = in[n * output_channel * in_h * in_w  + k * in_h * in_w + + (i + start_row_index) * in_w + j + start_col_index];
                }
            }
        }
    }
    
}

template<class I>
void matrixTrimPadding(I* in, I* out, size_t batch, size_t output_channel, size_t in_d, size_t in_h, size_t in_w, size_t out_d, size_t out_h, size_t out_w){
    if (in_h < out_h || in_w < out_w|| in_d < out_d){
        std::cout << "param error" << std::endl;
        std::cout << in_h << std::endl;
        std::cout << out_h << std::endl;
        std::cout << in_w << std::endl;
        std::cout << out_w << std::endl;
        std::cout << in_d << std::endl;
        std::cout << out_d << std::endl;
        throw "Error params in matrixTrimPadding! Please check!";
    }

    for (int n = 0 ; n < batch ; n++){
        for (int k = 0 ; k < output_channel ; k++){
            for (int d = 0 ; d < out_d ; d++){
                for (int i = 0 ; i < out_h ; i ++){
                    for (int j = 0 ; j < out_w ; j ++){
                        out[n * output_channel * out_d * out_h * out_w + k * out_d * out_h * out_w + d * out_h * out_w + i * out_w + j] 
                        = in[n * output_channel * in_d * in_h * in_w  + k * in_d * in_h * in_w + d * in_h * in_w + i* in_w + j];
                    }
                }
            }
        }
    }
}

/*
   在反向传播计算梯度后 进行裁剪

   可以推导出 一定有关系 ( outsize - insize ) % 2 = 0
*/

template<class I>
void matrixTrimBackData(I* in, I* out, size_t batch, size_t output_channel, size_t in_d, size_t in_h, size_t in_w, size_t out_d, size_t out_h, size_t out_w){
    int start_det_index = (in_d - out_d) / 2;
    int start_row_index = (in_h - out_h) / 2;
    int start_col_index = (in_w - out_w) / 2;

    for (int n = 0 ; n < batch ; n++){
        for (int k = 0 ; k < output_channel ; k++){
            for (int d = 0 ; d < out_d ; d++){
                for (int i = 0 ; i < out_h ; i ++){
                    for (int j = 0 ; j < out_w ; j ++){
                        out[n * output_channel * out_d * out_h * out_w + k * out_d * out_h * out_w + d * out_h * out_w + i * out_w + j] 
                        = in[n * output_channel * in_d * in_h * in_w  + k * in_d * in_h * in_w + (d + start_det_index) * in_h * in_w + (i + start_row_index)* in_w + j + start_col_index];
                    }
                }
            }
        }
    }
}

void matrixTrimRightWidth(Comp* in, Comp* out, size_t batch, size_t output_channel, size_t in_d, size_t in_h, size_t in_w, size_t out_w){
    for (int n = 0 ; n < batch ; n++){
        for (int k = 0 ; k < output_channel ; k++){
            for (int d = 0 ; d < in_d ; d++){
                for (int i = 0 ; i < in_h ; i ++){
                    for (int j = 0 ; j < out_w ; j ++){
                        out[n * output_channel * in_d * in_h * out_w + k * in_d * in_h * out_w + d * in_h * out_w + i * out_w + j] 
                        = in[n * output_channel * in_d * in_h * in_w  + k * in_d * in_h * in_w + d * in_h * in_w + i* in_w + j];
                    }
                }
            }
        }
    }
    std::cout << "matrixTrimRightWidth end" << std::endl;
    
}



// 向量叉乘，结果是一个向量
void vectorCrossMultip(Comp* a,Comp* b,Comp* c,size_t L) {
    for (int i = 0 ; i < L ; i++){
        c[i].real(a[i].real() * b[i].real() - a[i].imag() * b[i].imag());
        c[i].imag(a[i].imag() * b[i].real() + a[i].real() * b[i].imag());
    }
}


// 矩阵相加，将多个矩阵对应的位置相加，结果集中在指针的第一块位置
/**
    in-place transform
    batchN : 矩阵数量
*/
void matrixAddition(Comp* matrix, size_t batchN, size_t h, size_t w){
    for (int i = 1 ; i < batchN; i ++){
        for (int j = 0 ; j < h * w ; j++){
            matrix[j] += matrix[i * h * w +  j];
        }
    }
}

void matrixAddition(Comp* matrix, size_t batchN, size_t d, size_t h, size_t w){
    for (int i = 1 ; i < batchN; i ++){
        for (int j = 0 ; j < d * h * w ; j++){
            matrix[j] += matrix[i * d * h * w +  j];
        }
    }
}

void eliminateInputChannel(Comp* in , Comp* out , size_t h , size_t w){
    for (int i = 0 ; i < h * w ; i++){
       out[i] = in[i];
    }
}

void eliminateInputChannel(Comp* in , Comp* out , size_t d , size_t h , size_t w){
    for (int i = 0 ; i < d * h * w ; i++){
       out[i] = in[i];
    }
}

void oppositeImagPart(Comp* in, size_t length){
    for (int i = 0 ; i < length ; i++){
        in[i].imag(-in[i].imag()) ;
     }
}

void inputPadding(Real* in,Comp* out,size_t d, size_t h, size_t w , size_t p_d, size_t p_h, size_t p_w){

    for (int i = 0 ; i < d ; i++){
        for (int j = 0 ; j < h ; j++){
            for (int k = 0 ; k < w ; k++){
                out[(i + p_d) * (h + 2 * p_h) * (w + 2 * p_w) + (j + p_h) * (w + 2 * p_w) + (k + p_w)] 
                = in[i * h * w + j * w + k];
            }
        }
    }
}

void inputPadding(Real* in,Comp* out, size_t in_n, size_t in_c, size_t d, size_t h, size_t w , size_t p_d, size_t p_h, size_t p_w){

    for (int n = 0 ; n < in_n; n++){
        for (int m = 0 ; m < in_c; m++){
            for (int i = 0 ; i < d ; i++){
                for (int j = 0 ; j < h ; j++){
                    for (int k = 0 ; k < w ; k++){
                        out[n * in_c * (d + 2 * p_d) * (h + 2 * p_h) * (w + 2 * p_w) + m * (d + 2 * p_d) * (h + 2 * p_h) * (w + 2 * p_w) + (i + p_d) * (h + 2 * p_h) * (w + 2 * p_w) + (j + p_h) * (w + 2 * p_w) + (k + p_w)] 
                        = in[n * in_c * d * h * w + m * d * h * w +  i * h * w + j * w + k];
                    }
                }
            }
        }
    }
}

void inputPadding(Comp* in,Comp* out,size_t d, size_t h, size_t w , size_t p_d, size_t p_h, size_t p_w){

    for (int i = 0 ; i < d ; i++){
        for (int j = 0 ; j < h ; j++){
            for (int k = 0 ; k < w ; k++){
                out[(i + p_d) * (h + 2 * p_h) * (w + 2 * p_w) + (j + p_h) * (w + 2 * p_w) + (k + p_w)] 
                = in[i * h * w + j * w + k];
            }
        }
    }
}

// Because FFT computes a *one-sided* FFT, we need the final dimension to
// have *even* length.  Just pad with one more zero if the final dimension is odd.
void inputPaddingRightSide(Comp* in,Comp* out,size_t d, size_t h, size_t w , size_t p_d, size_t p_h, size_t p_w){

    for (int i = 0 ; i < d ; i++){
        for (int j = 0 ; j < h ; j++){
            for (int k = 0 ; k < w ; k++){
                out[i * h * (w + p_w) + j * (w + p_w) + k] 
                = in[i * h * w + j * w + k];
            }
        }
    }
}

int get_dim(size_t D,size_t H ,size_t W){
    if (W == 1)
    {
        return 0;
    } else if (D == H == 1)
    {
        return 1;
    } else if (D == 1)
    {
        return 2;
    } else {
        return 3;
    }
}

/*
    This method aim to rot a matrix 180 degree

    for example :

    before:

    [[1, 2, 3]
     [4, 5, 6]
     [7, 8, 9]]

    after:

    [[9, 8, 7]
     [6, 5, 4]
     [3, 2, 1]]
*/
void rot180Degree(Real* in, Real* out, size_t height, size_t width){
    for (int i = 0 ; i < height * width ; i++ ){
        out[i] = in[height * width - i - 1];
    }

}

void rot180Degree(Real* in, int length){
    size_t i = 0;
    size_t j = length - 1;
    Real tmp = 0;
    while (i < j)
    {
        tmp = in[i];
        in[i] = in[j];
        in[j] = tmp;
        i++;
        j--;
    }
}


void rot180Degree(Real* in, size_t n, size_t c, size_t d, size_t h, size_t w){
    for (int i = 0 ; i < n ; i++){
        for (int j = 0 ; j < c ; j++){
            rot180Degree(in + i * c * d * h * w + j * d * h * w , d * h * w);
        }
    }
}


void transposeByBlock(Real* in,Real* out,size_t n,size_t c,size_t h, size_t w){

    // for (int l = 0 ; l < n ;l++){
    //     for (int i = 0 ; i < c ; i++){
    //         for (int d = 0 ; d < h ; d++){
    //             for (int j = 0 ; j < w ; j++){
    //                   std::cout << in[l * c * h * w + i * h * w + d * w + j ] << "," ;
    //               }
    //           std::cout << std::endl;
    //         }
    //     std::cout  << std::endl; 
    //   }
    //   std::cout  << std::endl; 
    // }
    // std::cout  << std::endl;


    for (int i = 0 ; i < n ; i ++){
        for (int j = 0 ; j < c ; j ++){
            for (int k = 0 ; k < h ; k ++){
                for (int l = 0 ; l < w ; l ++){
                    out[j * n * h * w + i * h * w + k * w + l ] 
                    = in[i * c * h * w + j * h * w + k * w + l ];
                }
            }
        }
    }



    // for (int l = 0 ; l < n ;l++){
    //     for (int i = 0 ; i < c ; i++){
    //         for (int d = 0 ; d < h ; d++){
    //             for (int j = 0 ; j < w ; j++){
    //                   std::cout << out[l * c * h * w + i * h * w + d * w + j ] << "," ;
    //               }
    //           std::cout << std::endl;
    //         }
    //     std::cout  << std::endl; 
    //   }
    //   std::cout  << std::endl; 
    // }
    // std::cout  << std::endl;

}


void transposeByBlock(Real* in, Real* out, size_t n, size_t c, size_t d, size_t h, size_t w){

    // Real* tmp = (Real*) malloc(n * c * d * h * w * sizeof(Real));
    for (int i = 0 ; i < n ; i ++){
        for (int j = 0 ; j < c ; j ++){
            for (int m = 0 ; m < d ; m ++){
                for (int k = 0 ; k < h ; k ++){
                    for (int l = 0 ; l < w ; l ++){
                        out[j * n * d * h * w + i * d * h * w + m * h * w + k * w + l ] 
                        = in[i * c * d * h * w + j * d * h * w + m * h * w + k * w + l ];
                    }
                }
            }
        }
    }

    std::cout << "transposeByBlock end " << std::endl;
    // in = tmp;
    // free(tmp);

}
void transposeByBlockAndRot180(Real* in, Real* out, size_t n, size_t c, size_t d, size_t h, size_t w){

    for (int i = 0 ; i < n ; i ++){
        for (int j = 0 ; j < c ; j ++){
            for (int m = 0 ; m < d ; m ++){
                for (int k = 0 ; k < h ; k ++){
                    for (int l = 0 ; l < w ; l ++){
                        out[j * c * d * h * w + i * d * h * w + m * h * w + k * w + l ] 
                        = in[i * c * d * h * w + j * d * h * w + m * h * w + k * w + l ];
                    }
                }
            }
        }
    }

}
void copy_through_channel(Real* in, Real* out, size_t n, size_t c, size_t d, size_t h, size_t w){
        for (int i = 0 ; i < c ; i++ ){
            for (int j = 0 ; j < d ; j++ ){
                for (int k = 0 ; k < h ; k++ ){
                    for (int l = 0 ; l < w ; l++ ){
                        out[i*d*h*w + j*h*w + k*w + l] = in[j*h*w + k*w + l];
                    }
                }
            }
    }
}


// trim from start (in place)
static inline void ltrim(std::string &s) {
    s.erase(s.begin(), std::find_if(s.begin(), s.end(), [](int ch) {
        return !std::isspace(ch);
    }));
}

// trim from end (in place)
static inline void rtrim(std::string &s) {
    s.erase(std::find_if(s.rbegin(), s.rend(), [](int ch) {
        return !std::isspace(ch);
    }).base(), s.end());
}

// trim from both ends (in place)
static inline void trim(std::string &s) {
    ltrim(s);
    rtrim(s);
}




// // Non-square matrix transpose of matrix of size r x c and base address A
// void MatrixInplaceTranspose(int *A, int r, int c)
// {
//     int size = r*c - 1;
//     int t; // holds element to be replaced, eventually becomes next element to move
//     int next; // location of 't' to be moved
//     int cycleBegin; // holds start of cycle
//     int i; // iterator
//     bitset<HASH_SIZE> b; // hash to mark moved elements
  
//     b.reset();
//     b[0] = b[size] = 1;
//     i = 1; // Note that A[0] and A[size-1] won't move
//     while (i < size)
//     {
//         cycleBegin = i;
//         t = A[i];
//         do
//         {
//             // Input matrix [r x c]
//             // Output matrix
//             // i_new = (i*r)%(N-1)
//             next = (i*r)%size;
//             swap(A[next], t);
//             b[i] = 1;
//             i = next;
//         }
//         while (i != cycleBegin);
  
//         // Get Next Move (what about querying random location?)
//         for (i = 1; i < size && b[i]; i++)
//             ;
//         cout << endl;
//     }
// }

#endif
